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>(-lim),ifss<lim>(-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>(-lim), $
ifss<lim>(-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 ;=======================================================================