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,/ed),-width/2,-height/2)
Asmooc=smooth(abs(A),ns,/ed)
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 ;=======================================================================