pro trace_clean,data,cleaned_data,thresh
;+
; Project     : TRACE 
;
; Name        : TRACE_CLEAN
;
; Purpose     : Despiking of TRACE images
;
; Category    : Image processing
;
; Explanation : This algorithm cleans up spiky pixel in TRACE images which are either
;		high because of a cosmic ray hit or represent "hot pixels" that are 
;		temporarily enhanced. Because this algorithm is based on spatial and
;		temporal next neighbors (from the 3x3x3 cube around the spiky pixel) it 
;		requires a 3D data cube with a sequence of images DATA(*,*,NTIMES), 
;		NTIMES>3. It filters out noisy spikes that exceed a threshold factor THRES 
;		above the average of 2x8 nearest spatial neighbors of the preceding and
;		following image. The 8 nearest spatial neighbors in the simultaneous image
;		are not used because of the ring-like sidelobes produced by the jpeg compression
;		around a cosmic ray hit. The spiky pixel is replaced by the average of the
;		nearest 2 temporal pixels if the spikly pixel is a cosmic ray hit but 
;		not a "hot pixel". The noisy pixel is replaced by average of the 2x8 nearest 
;		(non-simultaneous) spatial neighbors if it is a "hot pixel".
;		This algorighm works similar to TRACE_DESPIKE and TRACE_DESTREAK, 
;		but seems to produce cleaner images for long exposure times where 
;		cosmic ray hit rate is high. Multiple iterations of this algorithm are 
;		recommended for deep cleaning (e.g. THESH=1.15), while a single iteration
;		if sufficient for coarse cleaning (e.g. THRESH=1.5).
;		This IDL algorithm is a vectorized version. 
;
; Syntax      : IDL> trace_clean,data,datac,thresh
;
; Examples    : IDL> read_trace,file,-1,index,data
;		IDL> thresh	=1.5
;               IDL> trace_clean,data,datac,thresh
;		IDL> tv,bytscl(data-datac)		;display difference	
;
; Inputs      : data	-  sequence of images [nx,ny] in form of data cube [nx,ny,nt], nt>3
;		thresh  -  threshold factor (e.g. 1.15-1.5) 
;
; Opt. Inputs : 
;
; Outputs     : datac	-  cleaned data cube
;               file='trace_clean.dat' containing statistics of cleaned pixels
;
; Opt. Outputs: None
;
; Keywords    : None 
;
; Common      : None
;
; Restrictions: None
;
; Side effects: None
;
; History     : Version 1,  30-MAR-1990,  Markus J. Aschwanden, LMSAL,  Written
;
; Contact     : aschwanden@lmsal.com
;-

if n_params(1) le 2 then thresh=1.5
print,'Threshold= ',thresh
cleaned_data=data
dim	=size(data)
ndim	=dim(0)
nx	=dim(1)
ny	=dim(2)
if (ndim eq 3) then nz=dim(3)
if (ndim eq 2) then nz=1
wdef,im=data(*,*,0)
openw,2,'trace_clean.dat',/append
printf,2,'NEXT ITERATION________________________________________________'
for iz=0,nz-1 do begin
 nclean	=0
 print   ,'Filter image #',iz,'   size=',nx,ny
 iz1	=iz-1	&if (iz eq 0)    then iz1=iz+2
 iz2	=iz+1	&if (iz eq nz-1) then iz2=nz-3
 for j=0,ny-1 do begin
  j1	=j-1	&if (j eq 0)    then j1=j+2
  j2	=j+1	&if (j eq ny-1) then j2=ny-3
  a1	=shift(data(*,j1,iz1),-1)	;8 nearest neighbors at t(i-1) 
  a2	=      data(*,j1,iz1)    
  a3	=shift(data(*,j1,iz1),+1) 
  a4	=shift(data(*,j ,iz1),-1) 
  a0	=      data(*,j ,iz1)    
  a5	=shift(data(*,j ,iz1),+1) 
  a6	=shift(data(*,j2,iz1),-1) 
  a7	=      data(*,j2,iz1)    
  a8	=shift(data(*,j2,iz1),+1) 
  b1	=shift(data(*,j1,iz2),-1) 	;8 nearest neighbors at t(i+1)
  b2	=      data(*,j1,iz2)    
  b3	=shift(data(*,j1,iz2),+1) 
  b4	=shift(data(*,j ,iz2),-1) 
  b0	=      data(*,j ,iz2)    
  b5	=shift(data(*,j ,iz2),+1) 
  b6	=shift(data(*,j2,iz2),-1) 
  b7	=      data(*,j2,iz2)    
  b8	=shift(data(*,j2,iz2),+1) 
  a_avg	=(a0+a1+a2+a3+a4+a5+a6+a7+a8)/8. ;average at time t(i-1)
  b_avg	=(b0+b1+b2+b3+b4+b5+b6+b7+b8)/8. ;average at time t(i+1)
  zav8	=((a_avg+b_avg)/2.)>0		 ;average of 16 pixels
  zav1	=((a0+b0)/2.)>0			 ;average of 2 nearest pixels in time
  zgood	=zav1
  ibad	=where(zav1 gt zav8*thresh,nbad)     ;bad high pixel
  if (nbad ge 1) then zgood(ibad)=zav8(ibad) ;replace bad 2-pixel average by 16-pixel average
  ispike=where(data(*,j,iz) gt zgood*thresh,nspike)	  ;spike detection
  if (nspike gt 0) then cleaned_data(ispike,j,iz)=zgood(ispike) ;replace central noispy spike
  nclean=nclean+nspike
 endfor
 print ,  'Fraction of noisy pixels cleaned = ',float(nclean)/float(nx*ny)
 printf,2,'image #',iz,'   Fraction of noisy pixels cleaned = ',float(nclean)/float(nx*ny)
 if (iz eq 0) then begin
  loadct,3
  statistic,cleaned_data(*,*,iz),zav,zsig
  z1	=zav-2*zsig
  z2	=zav+2.*zsig
 endif
 tv,bytscl(cleaned_data(*,*,iz),min=z1,max=z2)
endfor
close,2
end