FUNCTION star_squasher, fimage, ik, i0, nsigma=nsigma, mask_size=mask_size, $ debug=debug, edge_pix=edge_pix, nit=nit ; ========================================================================= ;+ ; PROJECT: ; Solar-B / XRT ; ; NAME: ; STAR_SQUASHER ; ; CATEGORY: ; Data calibration ; ; PURPOSE: ; This routine determines the locations of spikes, stars ; smudges & (remnant) streaks in the FFT of an image and suppresses ; (squashes!) these naughty features. The affected pixels are ; constrained to be <= the smoothed median value of nearby ; unaffected "background" pixels. ; Optionally outputs the pixel locations it has corrected (ik) ; and the pixel locations it has shielded from correction (i0) ; ; CALLING SEQUENCE: ; fft_out = STAR_SQUASHER (fft_in [,ik] [,i0]) ; ; INPUTS: ; FFT_IN - [Mandatory] (2-dim complex array, [Nx,Ny]) ; The FFT of an image containing read-out signals ; ; OUTPUTS: ; FFT_OUT - [Mandatory] (2-dim complex array, [Nx,Ny]) ; The FFT of an image with naughty read-out signals ; (mostly) removed ; IK - [Optional] (1-dim long array) ; Indicies of pixels in FFT_IN which were corrected ; I0 - [Optional] (1-dim long array) ; Indicies of pixels in FFT_IN which shielded from ; correction (contain significant data Fourier power) ; Note that some indicies appear in both I0 and IK; ; these were flagged as potential read-out signals, ; but most of their Fourier power was determined to ; arise from data ; ; ; EXAMPLES: ; ; Basic usage: ; IDL> fft_out = star_squasher(fft_in) ; ; You want the the indicies of changed/shielded pixels: ; IDL> fft_out = star_squasher(image_in, ik, i0) ; ; ; COMMON BLOCKS: ; None. ; ; KEYWORDS: ; NSIGMA - number of standard deviations ; beyond which an FFT pixel is ; considered bad. Default is 4.5 ; (Note: much below 4 not recommended) ; MASK_SIZE - size of neighborhood around ; k=0 which is not filtered. ; Default is 20 for a 512x512 image. ; EDGE_PIX - width in pixels of a border around ; the FFT that is not to be modified. The ; rationale is that typically, there is ; considerable power in the edges of the ; FFT of a rectangular image. Default = 5 ; for a 512x512 image. ; DEBUG - if set, give verbose diagnostic output. ; NIT - # of iterations; number of times to filter the FFT. ; It is MUCH faster and more precise to use ; NIT than to run star_squasher several ; times. caution - NOT TESTED. ; NOTES: ; 1) The amplitude spectrum abs(fft(image)) is broken ; down into rings of nearly constant k. The ; population in each ring is analyzed, and pixels ; that are more than NSIGMA standard deviations away from ; the local median/smoothed average are considered "bad". ; Corresponding pixels are replaced by (hopefully) appropriate ; background values. For finer points, see ; the list of keywords above. ; 2) For images of sizes 128xN or smaller, some of the Fourier ; signals may be incompletely removed. ; ;MODIFICATION HISTORY: progver = 'v2007.Apr.23' ; --- (SS, JC, MW) written (as xrt_rnoise) ; by JC (1/2007), adapted from ; trace_rnoise.pro. Considerably modified ; by SS (4/2007). Cleaned up by MW. progver = 'v2007.Nov.27' ; --- (SS) Fixed bug for small images (side=128). ; ; ;- ; ========================================================================= isize=size(fimage) width=isize(1) height=isize(2) if not keyword_set(nsigma) then nsigma=4.5 if not keyword_set(mask_size) then mask_size=20*max([width,height])/512. if not keyword_set(nit) then nit=1 f=fimage ;===== nominal smoothing size for 512x512 image nsmoo=11. ;===== minimum image size for non-tested ring searches minsize=128 ;===== make k matrix of distances from 0,0 k = dist(width,height) ;===== calculate the logarithm of k scaled so that k changes by 1 per pixel ;===== in the neighborhood of k=msize. Truncate to integer. kmap = fix(alog(k/float(mask_size))/alog(1.0+1.0/mask_size)) ;===== prepare for edge exclusion via keyword edge_pix if (n_elements(edge_pix) eq 0) then edge_pix=5 if (edge_pix gt 0) then begin kmap(0:edge_pix-1, *)=0 kmap(width-edge_pix:width-1, *)=0 kmap(*, 0:edge_pix-1)=0 kmap(*, height-edge_pix:height-1)=0 endif ;===== save indicies which are protected from correction i0=where(kmap le 0) N = max(kmap) if keyword_set(debug) then print,'constructed kmap with ',N,' levels.' nbad=0L ;initialize bad pixel count if keyword_set(debug) then f_old=fimage ;===== Fourier amplitudes A = abs(f) ;===== scale smoothing size ns=round(nsmoo*max([width,height])/512) ;===== smoothed transform (smooth edges also) Asmoo=shift(smooth(abs(shift(A,width/2,height/2)),ns,/edge_truncate),-width/2,-height/2) Asmooc=smooth(abs(A),ns,/edge_truncate) Asmoo(width/2-ns:width/2+ns,*)=Asmooc(width/2-ns:width/2+ns,*) Asmoo(*,height/2-ns:height/2+ns)=Asmooc(*,height/2-ns:height/2+ns) mn=min(Asmoo) csig=-0.5 ; -1 ;===== setup index vector of corrected points ik=-1 for iteration=1,nit do begin ;===== check image size - if length & width > minsize do normal ring search if width gt minsize and height gt minsize then begin for k = 1, N-1 do begin ;N-1 in case too few pixels with k=N. ;===== for given ring kmap=k do statistics ss = where(kmap eq k) mom = moment(A(ss)) A_med = median(A(ss)) A_sigma = sqrt(mom(1)) ;===== find where Fourier amp > median_ring + nsigma*sigma_ring sss = where(A(ss) gt A_med + nsigma*A_sigma) if sss(0) ne -1 then begin ss = ss(sss) ik=[ik,ss] rfss=float(f(ss)) ifss=imaginary(f(ss)) ;===== set limit vector to be the lesser of median_ring or a local smooth ;===== plus csig*sigma_ring; insist this limit is > smoothed global min lim=(min([[Asmoo(ss)],[ss*0.+A_med]],dim=2) + $ csig*A_sigma)>mn ;==== constrain real and imaginary parts to be < limit, retaining sign f(ss)=complex(rfss(-lim),ifss(-lim)) nbad = nbad + n_elements(ss) endif endfor ;===== if length or width <= minsize do ring search but check number of ring ;===== elements endif else begin for k = 1, N-1 do begin ;N-1 in case too few pixels with k=N. ;===== for given ring kmap=k do statistics ss = where(kmap eq k,nss) if nss ge 2 then begin mom = moment(A(ss)) A_med = median(A(ss)) A_sigma = sqrt(mom(1)) ;===== find where Fourier amp > median_ring + nsigma*sigma_ring sss = where(A(ss) gt A_med + nsigma*A_sigma) if sss(0) ne -1 then begin ss = ss(sss) ik=[ik,ss] rfss=float(f(ss)) ifss=imaginary(f(ss)) ;===== set limit vector to be the lesser of median_ring or a local smooth ;===== plus csig*sigma_ring; insist this limit is > smoothed global min lim=(min([[Asmoo(ss)],[ss*0.+A_med]], $ dim=2) + csig*A_sigma)>mn ;==== constrain real and imaginary parts to be < limit, retaining sign f(ss)=complex(rfss(-lim), $ ifss(-lim)) nbad = nbad + n_elements(ss) endif endif endfor end endfor if n_elements(ik) gt 1 then ik=ik(1:*) if keyword_set(debug) then begin print,'suppressed ',nbad,' "bad" FFT pixels.' xtv, float(f ne f_old) print,'Display shows FFT bad pixel map.' stop endif result=f return, result END ;=======================================================================