;+ ;NAME: SOT_NOSPIKE.PRO ; ;PURPOSE: Despike an image using a median filter / dilation ; This is a modified version of "nospike.pro" for SOT data ; ;Method: spikes wtih ; * Brighter than (`threshold' percent AND average) ; OR ; * Darkter than ('threshold' percent AND average) ; Fills hits with median of neighbouring pixels that were not ; flagged as radiation hits. ; Contrary to trace_despike, this routine often doesn't need ; iterations, and two iterations usually suffices. ;SUBROUTINES: ; wv_a_trous ; ;CALLING SEQUENCE: ; image_out = sot_nospike(index, image) ; or image_out = sot_nospike(image) ; ;OUTPUT: returns despiked image array ; ;CALLS: ; avg, label_region ; ;OPTIONAL KEYWORD INPUT: ; ; nspikes - number of spikes (not pixels) removed ; imap - map of corrected pixels ; threshold - fractional threshold for brightness (def: 0.15) ; brightthreshold- fractional threshold for brightness in bright regions ; (def: 0.15) ; minimum - do not correct BRIGHT pixels dimmer than minimum ; (def: avg(image)) ; maximum - do not correct DARK pixels brighter than maximum ; (def: avg(image)) ; flare - sets brightthreshold to 0.3 ; silent - doesn't print out messages ; update - If set, update history record in INDEX ; ;Restrictions: ; the call to RSI 'dilate' varies slightly between pre/after V5.3 ; (RSI added '/constrain' and 'background' keywords in 5.3) ; ;HISTORY: ; ; Bart De Pontieu (BDP) 27 October 2000. ; 31-October 2000, S.L.Freeland rename Bart's ; 'trace_nospike.pro' to 'nospike.pro' and add it ; to the growing SSW/gen CCD cleanup arsenal ; added version dependent 'dilate' call ; Zarro (EER/GSFC) 28 May 2003 - added "where" check ; Y. Katsukawa 16 Apr 2007 Ported for despike of HINODE SOT data ; Dark pixel correction is added. ; ;v1.1 Adapted for use with Hinode/SOT NFI images. T. Berger 25-Oct-2007. ;- PRO wv_a_trous, cin,cout,wout,j_max ;; Fast version of algorithm described on p.24 and following ;; in "Image Processing and Data Analysis" by Starck et al. ;; Similar to unsharp masking scin = size(cin) nx = scin(1) ny = scin(2) cin2 = float(cin) cout = fltarr(nx,ny,j_max+1) cout(0,0,0) = cin2 wout = cout*0. FOR j=1,j_max DO BEGIN n = 2^(j-1) kernel = FLTARR(2*n+1,2*n+1)*0. kernel[n*INDGEN(3),0] = [0.0625,0.125,0.0625] kernel[n*INDGEN(3),n] = [0.125,0.25,0.125] kernel[n*INDGEN(3),2*n] = [0.0625,0.125,0.0625] cout[0,0,j] = CONVOL(cout[*,*,j-1],kernel,/edge_truncate) wout[0,0,j] = cout[*,*,j-1]-cout[*,*,j] ENDFOR return END ;----------------------------------------------------------------------------- PRO sot_nospike, index, data, index_out, data_out, $ nspikes=nspikes, imap=imap, $ threshold=threshold, brightthreshold=brightthreshold, $ minimum=minimum, maximum=maximum, $ version=progver, name=prognam, $ loud=loud, verbose=verbose, $ flare=flare, update=update ;----------------------------------------------------------------------------- progver = 'V1.1' prognam = 'SOT_NOSPIKE.PRO' time0 = SYSTIME(1) retrn = 0 loud = KEYWORD_SET(loud) verbose = KEYWORD_SET(verbose) if (verbose eq 1) then loud = 1 if (loud) then print, 'Start ',prognam,' ',progver if N_PARAMS() lt 2 then begin MESSAGE,/INFO,'Useage: data_out = sot_nospike(index,data,/KEYWORDS...)' data_out = data RETURN END ; replacement thresholds IF KEYWORD_SET(threshold) THEN threshold=threshold ELSE threshold=0.15 IF KEYWORD_SET(flare) THEN threshold2 = 0.3 IF KEYWORD_SET(brightthreshold) THEN threshold2 = brightthreshold ELSE BEGIN IF NOT(keyword_set(flare)) THEN threshold2 = 0.15 ENDELSE ds = SIZE(data) nx = ds[1] ny = ds[2] case ds[0] of 1: begin MESSAGE,'2- or 3-D image required. Returning...' data_out = data RETURN end 2: nimage = 1 3: nimage = ds[3] else: begin MESSAGE,'2- or 3-D image arrays only. Returning...' data_out = data RETURN end end ; Dilation Kernel ker = INTARR(3,3) ker(0,1) = 1 ker(1,*) = 1 ker(2,1) = 1 data_out = data index_out = index for i=0,nimage-1 do begin image = data[*,*,i] image_type = fg_image_type(index,i) tp = data_type(image) ; minimum intensity of radiation hit IF KEYWORD_SET(minimum) THEN minimum=(minimum > 0) ELSE minimum=avg(image) ; maximum intensity of pixel defects IF KEYWORD_SET(maximum) THEN maximum=(maximum > 0) ELSE maximum=avg(image) switch image_type of 'FG': 'SI': begin ; periodic extension to deal with hits at the edge image2 = MAKE_ARRAY(nx+4,ny+4,type=tp) image2(2,2) = image image2(0,0) = image2(4,*) image2(1,0) = image2(3,*) image2(nx+3,0) = image2(nx-1,*) image2(nx+2,0) = image2(nx,*) image2(0,0) = image2(*,4) image2(0,1) = image2(*,3) image2(0,ny+3) = image2(*,ny-1) image2(0,ny+2) = image2(*,ny) th = FLTARR(nx+4,ny+4)+threshold ; Call `a trous algorithm for unsharp masking wv_a_trous,image2,cim,wim,1 ; Define bright mask ci1 = SMOOTH(image2,11,/edge_truncate) aci1 = avg(ci1) sci1 = stdev(ci1) ;-- changed by DMZ chk=WHERE(ABS(ci1-aci1) GE sci1, count) if count gt 0 then th(chk)=threshold2 ; Define minimum brightness mask wi1 = REFORM(wim(*,*,1)) wi1b = wi1*0. wmin = WHERE(image2 GE minimum, nmin1) IF nmin1 GT 0 THEN BEGIN wi1b(wmin) = wi1(wmin)/FLOAT(image2(wmin)) ENDIF ; Define maximum brightness mask wi1a = wi1*0. wmax = WHERE(image2 LE maximum, nmax1) IF nmin1 GT 0 THEN BEGIN wi1a(wmax) = wi1(wmax)/(FLOAT(image2(wmax))>10) ENDIF ; Dilate & threshold unsharp mask mask1 = ((wi1b - th) GE 0) or ((wi1a + th) LE 0) skl = DILATE(mask1,ker) ; Label all radiation hits in dilated mask b_curv = LABEL_region(skl,/eight) h_curv = HISTOGRAM(b_curv,reverse_indices=r) estring=(['b_curv2 = dilate(b_curv,ker,/gray)', $ 'b_curv2 = dilate(b_curv,ker,/gray,/ulong,/constrained)']) estat=execute( (estring)(since_version('5.3'))) estring=(['b_curv3=dilate(b_curv2,ker,/gray)',$ 'b_curv3 = dilate(b_curv2,ker,/gray,/ulong,/constrained)']) estat=execute( (estring)(since_version('5.3'))) h_curv3 = histogram(b_curv3-b_curv,reverse_indices=r3) curv_tot2 = skl*0 skl_tot = skl*0. ; Fill in labeled radiation hits FOR kk=1L,n_elements(h_curv)-1L DO BEGIN p = r(r[kk]:r[kk+1]-1) p3 = r3(r3[kk]:r3[kk+1]-1) curv_tot2[p] = 1 skl_tot[p] = median(image2[p3]) ENDFOR ; Clean up nspikes = LONG(N_ELEMENTS(h_curv)) imap = curv_tot2(2:nx+1,2:ny+1) ; Messages IF (loud) THEN BEGIN MESSAGE, /info, 'Corrected '+STRCOMPRESS(STRING(N_ELEMENTS(h_curv)) $ ,/remove_all)+' radiation hits and bad pixels (>'+$ STRCOMPRESS(STRING(FIX(threshold*100)),/remove_all)+' %)' ENDIF ; Calculate new image skl_final = make_array(nx,ny,type=tp) skl_final[0,0] = (skl_tot*curv_tot2+image2*(1-curv_tot2))[2:nx+1,2:ny+1] image = skl_final BREAK end ;FG and SI case 'SQ': 'SU': 'SV': begin MESSAGE,/INFO,'Stokes Q, U, and V images not despiked...' BREAK end else: begin MESSAGE,/INFO,'Unknown image type. Returning...' data_out = data RETURN end endswitch data_out[*,*,i] = image endfor run_time = SYSTIME(1) - time0 if (loud) then PRINT, 'SOT_NOSPIKE took '+strtrim(run_time,2)+' seconds for 1 image' ; Update history record in the index if requested if KEYWORD_SET(update) and (N_PARAMS() ge 2) then begin tagval = 'Radiation despiking applied.' update_history,index_out, tagval, version=progver if (loud) then MESSAGE, /INFO, 'HISTORY record updated: '+tagval endif RETURN END