; $Id: $SSW/proba2/swap/idl/prep/p2sw_despike.pro#1 $ ; ; Copyright (c) Trinity College Dublin. All rights reserved. ;+ ; NAME: ; P2SW_DESPIKE ; ; PURPOSE: ; This function applies despiking of normal cosmic ray hits ; and locally dark pixels from PROBA2/SWAP images. ; ; CATEGORY: ; PROBA2/SWAP calibration. ; ; CALLING SEQUENCE: ; ; Result = P2SW_DESPIKE( Index, Image [, SPIKE_SIG=spike_sig, $ ; COR_MASK=cor_mask, VERSION=version, $ ; RUNID=runid, /VERBOSE, /LMAT, $ ; /APPLY_GENTLE_FILTER, _EXTRA=ex ] ) ; ; INPUTS: ; Index: FITS header in structure format for PROBA2/SWAP image ; to be processed. ; ; Image: PROBA2/SWAP image to be processed. ; ; OPTIONAL INPUTS: ; SPIKE_SIG: Number of local standard deviations above the local ; median to identify cosmic ray pixels. DEFAULT: 6. ; ; VERSION: Master P2SW_PREP code version number for LMAT logging ; and HISTORY tag entries. DEFAULT: null string (''). ; ; RUNID: Unique processing run identifier for P2SC pipeline ; operations. Used only if LMAT is set. DEFAULT: 0L. ; ; KEYWORD PARAMETERS: ; VERBOSE: Set this keyword to print run-status information to ; standard output ; ; LMAT: Set this keyword to enable default reporting mode for ; P2SC pipeline operations. Not for use via SSW. ; ; APPLY_GENTLE_FILTER: Set this keyword to apply a filter that ; restricts the value of every pixel to lie ; between the min and max values of its eight ; neighbors. This function is useful for ; filtering images with few, strong outlier ; pixels while preserving the overall ; sharpness of the image. ; ; OUTPUTS: ; This function returns the input image corrected for normal ; cosmic ray hits. Input Index structure is updated during ; processing - e.g., HISTORY appended. ; ; OPTIONAL OUTPUTS: ; COR_MASK: Output is passed as binary mask identifying corrected ; pixels (1) and uncorrected pixels (0). ; ; PROCEDURE: ; Input image has a 5x5 pixel median filter applied and each ; pixel has a local standard deviation calculated using the ; Median Absolute Deviation method within the 5x5 pixel locale ; (st. dev. ~ 1.483*MAD; REFERENCES). ; Cosmic ray affected pixels are identified as having values ; greater than SPIKE_SIG * st. dev. above the local median. Those ; pixels selected are replaced by their 5x5 pixel median-filter ; value. Locally dark pixels are also caught by this method, due ; to the testing of absolute deviation from the 5x5 median value. ; ; MODIFICATION HISTORY: ; Written by: D. Shaun Bloomfield, 24-Jan-2009. ; Apr 2009 Added IDL-format description. ; Jun 2009 Moved FITS history updating into sub-routines. ; Apr 2010 Changed to standard deviation calculation by MAD. ; Added COR_MASK output mask of corrected pixels. ; May 2010 Finalized IDL-format description. ; Sep 2010 Fixed problem that zeroed corners in subfields (dbs) ; Sep 2012 Fixed problem that skipped pixels on image edges (dbs) ; Jul 2016 Added apply_gentle_filter to allow for image filtering (edh) ; Jun 2017 Fixed code to calculate MAD on original image instead of ; on partially despiked image (dbs, edh) ; Jun 2017 Speed up processing by replacing FOR-loop with data cube (dbs, edh) ;- FUNCTION p2sw_despike, index, image, $ spike_sig=spike_sig, version=version, $ verbose=verbose, $ apply_gentle_filter = apply_gentle_filter, $ cor_mask=cor_mask, _extra=ex prognam = 'p2sw_despike.pro' ;--- Check keyword inputs and set keyword-dependencies IF KEYWORD_SET(spike_sig) THEN spike_sig = spike_sig ELSE spike_sig = 6. IF KEYWORD_SET(version) THEN version = version ELSE version = '' verbose = KEYWORD_SET(verbose) image1 = image subimage = image1(index.firstcol - 1: index.last_col - 1, index.firstrow - 1: index.last_row - 1) cor_mask = BYTE(image*0) ; extract subfield, use filter_image to ensure median catches spikes at ; the edge of the subfield as well xmax = (index.last_col - index.firstcol) ymax = (index.last_row - index.firstrow) med_im = FILTER_IMAGE( subimage, median = 5, /all_pixels ) dev_im = ABS( subimage - med_im ) sub_stack = abs([ [[shift(subimage, -2, -2) - med_im]], $ [[shift(subimage, -2, -1) - med_im]], $ [[shift(subimage, -2, 0) - med_im]], $ [[shift(subimage, -2, 1) - med_im]], $ [[shift(subimage, -2, 2) - med_im]], $ [[shift(subimage, -1, -2) - med_im]], $ [[shift(subimage, -1, -1) - med_im]], $ [[shift(subimage, -1, 0) - med_im]], $ [[shift(subimage, -1, 1) - med_im]], $ [[shift(subimage, -1, 2) - med_im]], $ [[shift(subimage, 0, -2) - med_im]], $ [[shift(subimage, 0, -1) - med_im]], $ [[shift(subimage, 0, 0) - med_im]], $ [[shift(subimage, 0, 1) - med_im]], $ [[shift(subimage, 0, 2) - med_im]], $ [[shift(subimage, 1, -2) - med_im]], $ [[shift(subimage, 1, -1) - med_im]], $ [[shift(subimage, 1, 0) - med_im]], $ [[shift(subimage, 1, 1) - med_im]], $ [[shift(subimage, 1, 2) - med_im]], $ [[shift(subimage, 2, -2) - med_im]], $ [[shift(subimage, 2, -1) - med_im]], $ [[shift(subimage, 2, 0) - med_im]], $ [[shift(subimage, 2, 1) - med_im]], $ [[shift(subimage, 2, 2) - med_im]] ]) mad_im = median(sub_stack, dim = 3) cor_mask = dev_im gt spike_sig * 1.483 * mad_im subimage[where(cor_mask)] = med_im[where(cor_mask)] ; Old method is inefficient ; subimage_original = subimage ; FOR x=0, xmax DO $ ; FOR y=0, ymax DO BEGIN ; local = subimage_original[ (x-2) > 0 : (x+2) < xmax, (y-2) > 0 : (y+2) < ymax ] ; MAD = MEDIAN( ABS( local - MEDIAN( local ) ) ) ; IF ( dev_im[x, y] GT (spike_sig*1.483*MAD) ) THEN BEGIN ; subimage[x, y] = med_im[x, y] ; cor_mask[x, y] = 1 ; ENDIF ; ENDFOR image1(index.firstcol - 1: index.last_col - 1, index.firstrow - 1: index.last_row - 1) = subimage ;--- Set HISTORY update text tagval = STRING( 'Despiked at the ', spike_sig, '-sigma level ('+STRTRIM( ROUND( TOTAL(cor_mask) ), 2 )+' pix)', FORMAT='(a,f3.1,a)' ) ;--- Report if requested IF (verbose) THEN MESSAGE, 'HISTORY record updated: '+tagval, /CONTINUE ;--- Add update text into HISTORY field of index structure UPDATE_HISTORY, index, tagval, CALLER='P2SW_DESPIKE '+version ; if extra filtering is requested, apply it here IF (apply_gentle_filter) THEN BEGIN image1 = gentle_filter(image1, /second_max) ;--- Set HISTORY update text tagval2 = STRING('Gentle filter applied') ;--- Report if requested IF verbose THEN MESSAGE, 'HISTORY record updated: '+tagval2, /CONTINUE ;--- Add update text into HISTORY field of index structure UPDATE_HISTORY, index, tagval2, CALLER='P2SW_DESPIKE '+version ENDIF RETURN, image1 END